% Created 2021, Mike Siedlik
% Free for research and non-commercial use
% Please cite Siedlik and Issadore, 2021

% This MATLAB code was created with the following version of MATLAB:
% MATLAB Version: 9.10.0.1739362 (R2021a) Update 5
% This code uses the following MATLAB toolboxes:
% Image Processing Toolbox

% -------------------------------------------------------------
% -------------------------------------------------------------
% SPECIFY DATA

% An example data set is provided as an image sequence, though the code
% could be easily modified to read videos or individual images.
data_name = 'A_ExampleData.tif';
% Read image sequence to obtain image dimensions and number of frames, and
% to preallocate a matrix to store image data
INFO = imfinfo(data_name);
% Filename to be used for the reconstructed, cropped image sequence
output_name = 'A_ExampleData-Reordered.tif';


% -------------------------------------------------------------
% -------------------------------------------------------------
% USER-DEFINED IMAGE PROCESSING PARAMETERS

% "Y-distance" above and below the channel. Since the example data uses
% a straight channel, this is used to identify the channel region for
% object detection
dy = 109;

% Threshold value for intensity-based thresholding and object detection
thresh = 145;

% Boundary region to use for identifying a droplet. In this application,
% the second entry is the most important, as the logic of this script is to
% identify the "last droplet within some bounds" as the reference droplet.
x = [225,525];

% Mean filter to be used for blurring images
h = 1/numel(ones(7))*ones(7);

% Box used to crop image (possibly unique to this example). For
% convenience, cropping is included here, where the frames and bounding box
% are determined by evaluating an image in the sequence.
% cB(1): top left corner of box to crop, x position
% cB(2): top left corner of box to crop, y position
% cB(3): box width
% cB(4): box height
cB = [266,57,278,139];
% Range of frames to store image data
frame_range = [275,875];



% -------------------------------------------------------------
% -------------------------------------------------------------
% -------------------------------------------------------------
% -------------------------------------------------------------
% IDENTIFY REFERENCE DROPLET IN EACH IMAGE

% Preallocate matrix to store reference droplet information
refData = zeros(numel(INFO),5);
% refData(:,1): frame number
% refData(:,2): object identifier (object in frame that is determined to be
% the reference droplet)
% refData(:,3): droplet area (pixels)
% refData(:,4): droplet x centroid position
% refData(:,5): droplet y centroid position


% Loop through each image. For each slice, the goal is to determine the x
% position of the centroid of the droplet between the x limits (column
% positions) in that image
for k = 1:numel(INFO)
    
    % 1) Load image frame
    im_k = imread(data_name,k);
    
    % 2) Consider only green channel of RGB image (arbitrary, could also
    % find mean by averaging channels).
    im_k = im_k(:,:,2);
    
    % 3) Apply mean filter to smooth intensity based information
    im_k = filter2(h,im_k);
    
    % 4) "Block out" "y" regions of image that do not correspond to the
    % channel
    im_k([1:dy,end-dy+1:end],:) = 255;
    
    % 5) Apply intensity based thresholding. Since this image data is based
    % on imaging droplets via transmitted light imaging that contain an
    % absorbing dye, create a binary image of all pixels LESS than the
    % specified thresshold.
    im_bw = im_k < thresh;
    
    % Temporary vector to store all detected droplets/objects
    temp_DropletInfo = [];
    
    % 6) Determine all objects in binarized image
    [L,num] = bwlabel(im_bw);
    
    % For each identified object, analyze the size to determine if it is a
    % droplet. % NOTE: other parameters, such as Circularity or AxisLength,
    % could be used as well. For intensity-based threshold of droplets
    % containing dye, Area is often sufficient.
    for m = 1:num
        
        % Obtain image corresponding to the m-th object and fill any holes
        im_obj = imfill(L == m,'holes');
        
        % Obtain information about the m-th reference object.
        objStats = regionprops(im_obj,'Area','Centroid');

        % Identify full droplets (as opposed to stray pixels or partial
        % droplets) based on object area.
        if objStats.Area > 2250
            
            % Store that droplet's information 
            temp_DropletInfo(end+1,1:5) = [k,m,objStats.Area,objStats.Centroid];

        end
    end
    
    % Store information for the "reference droplet," i.e. the "last droplet
    % within some bounds." Note: based on the MATLAB image identification
    % algorithm, objects IN THIS APPLICATION are by default ordered by x
    % position. If this was not the case, the temp_DropletInfo matrix could
    % be first sorted in the process of obtaining the last droplet within
    % the specified bounds.
    temp_DropletInfo = temp_DropletInfo(temp_DropletInfo(:,4) < x(2),:);
    refData(k,1:5) = temp_DropletInfo(end,:);

end

% -------------------------------------------------------------
% SAVE RESORTED FRAMES IN NEW IMAGE SEQUENCE OR VIDEO

% Sort reference droplet array data by x position of reference droplet
[~,I2] = sort(refData(:,4),'Ascend');
refData = refData(I2,:);


% Save each image frame in an image sequence. Note this saves a cropped
% image sequence (in XY and in T) based upon the parameters specified at
% the beginning.
for n = frame_range(1):frame_range(2)
    % Read the given frame
    im_n = imread(data_name,refData(n,1));
    
    % Write cropped frame
    imwrite(im_n(cB(2):cB(2)+cB(4),cB(1):cB(1)+cB(3),:),...
        output_name,'WriteMode','append')
end


